terraThis workshop will provide you with an introduction to manipulating
raster spatial data using the R package terra.
terra and its predecessor, raster are widely
used for spatial data manipulation in R. No prior experience of spatial
data is assumed, but this short workshop will not have time to delve
into some important aspects of spatial data such as projections.
terra’s predecessor raster had many of the
same functions as terra, and I will mention how functions
have changed or been renamed which might be helpful for people migrating
from using raster. In the words of the creator of both
packages: terra is simpler, faster and can do more, so
definitely switch to terra if you are still using
raster!
raster users note: Boxes like this will highlight
differences between the raster package and
terra.
Resources:
terra tutorial page is https://rspatial.org/spatial/index.htmlterra is
https://rspatial.github.io/terra/reference/terra-package.htmlterra and sf packages: https://mgimond.github.io/Spatial/index.htmlYou will need the terra package installed, which can be
done by running install.packages("terra"). If you have
problems, there are more details about installing it here. We can then
load it
library(terra)
#> terra 1.7.71
There are basically two types of spatial data: vector and raster
Can be points, lines or polygons. Useful for representing things like survey locations, rivers, and boundaries.
pts <- rbind(c(3.2,4), c(3,4.6), c(3.8,4.4), c(3.5,3.8), c(3.4,3.6), c(3.9,4.5)) |>
vect()
lnes <- as.lines(vect(rbind(c(3,4.6), c(3.2,4), c(3.5,3.8)))) |>
rbind(as.lines(vect(rbind(c(3.9, 4.5), c(3.8, 4.4), c(3.5,3.8), c(3.4,3.6)))))
lux <- vect(system.file("ex/lux.shp", package = "terra"))
par(mfrow = c(1,3))
plot(pts, axes = F, main = "Points")
plot(lnes, col = "blue", axes = F, main = "Lines")
plot(lux, "NAME_2", col = terrain.colors(12), las = 1, axes = F, main = "Polygons")
par(mfrow = c(1,3))
Raster data is a grid of rectangles, normally called cells. Each cell has a value, making rasters useful for storing continuous data, such as temperature and elevation.
Here is an example of raster data, where each cell in the raster represents elevation.
elev <- system.file("ex/elev.tif", package = "terra") |>
rast() |>
aggregate(fact = 2)
plot(elev, las = 1, main = "Elevation map")
elev |>
as.polygons(aggregate = FALSE, na.rm = FALSE) |>
lines(col = "grey40", lwd = 0.2)
Lets start by creating our own raster. We can create rasters from
scratch and load them from a file using the function
rast(). We can create a simple raster by specifying the x
and y limits for the raster and the resolution (how big each cell
is).
raster users note: rast() replaces
raster()
#create raster
ras <- rast(xmin = 0, xmax = 10, ymin = 0, ymax = 10, resolution = 2)
#see what we've created
ras
#> class : SpatRaster
#> dimensions : 5, 5, 1 (nrow, ncol, nlyr)
#> resolution : 2, 2 (x, y)
#> extent : 0, 10, 0, 10 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84
The figure below shows what most of the terms above refer to. As you can see, you don’t need to use all the terms to define a raster. Couple of other points:
data.frame and
as you can see, rasters in terra are of class
SpatRaster.rast() which coordinate reference
system to use, so it defaults to using longitude latitude coordinates,
also known as EPSG 4326.raster users note: rasters are now
SpatRaster class not RasterLayer
But what does the raster we created actually look like when plotted.
Lets see. All we need is plot()
plot(ras)
Why is there no plot? Because the raster we created is empty; there are no values associated with the the cells. Lets assign some values to each cell in the raster and try again. First we will find out how many cells are in our raster using `ncell()
ncell(ras)
#> [1] 25
Ok, now we know this lets give our raster cells values from 1 to 25:
values(ras) <- 1:25
plot(ras)
Now our raster has values, we get a plot! Each cell has an integer value between 1 and 25, with cell values increasing from left to right and top to bottom. So the values start being “filled up” in the top left, and finish in the bottom right.
Lets have another look at our raster properties
ras
#> class : SpatRaster
#> dimensions : 5, 5, 1 (nrow, ncol, nlyr)
#> resolution : 2, 2 (x, y)
#> extent : 0, 10, 0, 10 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84
#> source(s) : memory
#> name : lyr.1
#> min value : 1
#> max value : 25
We can now see a few extra pieces of information compared to last time:
sources(s): where is the data held on your computer? It
says memory for this raster, indicating that the raster is
in the computer memory. Rasters can also be held on your hard disk, in
which case this will be the file name of the raster. We won’t go into
details here, but terra is smart about loading data into
memory, only doing so when it needs to and it thinks it will have enough
space.name: what is the raster called?min value & max value: the minimum and
maximum values in the rasterOk, now we understand the basic structure of a raster, lets load some data.
Sea surface temperature can be measured by satellites, and the National Oceanic and Atmospheric Administration (NOAA) in the U.S. has been doing so since 1972! Daily sea surface temperature for the world from 1985 to the present is available via the NOAA Coral Reef Watch website.
We are going to explore sea surface temperature (SST) data for the
Great Barrier Reef in Australia. All the data is in the
data folder in the Github repository, and full details on
how I got the data are in the data_prep.R script in that
folder.
We can use the same rast() command that we used to make
our own raster to load data The data is stored in GeoTiff (.tif) format,
which is widely used for raster data. You can get a full list of data
formats that terra can read and write by running
gdal(drivers = TRUE).
sst <- rast("../data/gbr_temp_2023_05_31.tif") #load the data
sst #view the properties
#> class : SpatRaster
#> dimensions : 801, 1101, 1 (nrow, ncol, nlyr)
#> resolution : 0.04999999, 0.05 (x, y)
#> extent : 110, 165.05, -45, -4.95 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)
#> source : gbr_temp_2023_05_31.tif
#> name : CRW_SST_2
#> min value : 9.59
#> max value : 30.46
#> unit : Celsius
If the data loaded correctly, you should see the properties as shown above. We have all the properties I described earlier, but there are a couple of things worth noting:
source is the filename of the data we loaded. This
means it is on disk, not in memory. We can double check by running
inMemory(sst) which should return FALSE.unit is Celsius. Most rasters don’t come with the
units, so you will normally need to look at metadata to find this
information.We want to plot our data to see what it looks like. This is the first thing you should do with almost any spatial data to do a quick check that it looks right
plot(sst)
The values in the legend seem about what we would expect for water temperatures, ranging from ~10 - 30\(^\circ\) Celsius.
The white areas of the plot are cells with NA values,
which in this case are land; there are no sea surface temperature data
for the land! We can see the land mass of Australia taking up much of
the map and Papua New Guinea to the north.
We can set the NA values to a different colour if we want, which can be helpful if you want to see some NA cells that are getting lost against the other colours.
plot(sst, colNA = "black")
Note that every cell within our raster has to have a value.
To put our raster data in context, lets plot the Great Barrier Reef
Marine Park boundary (dowloaded from here).
This is vector data. You might be used to handling vector data with the
sf package, but terra can also be used for
vector data manipulation. We can load a vector using
vect()
raster users note: terra has its own
methods for handling vector data unlike raster which used
the sp package for vector data handling. Vector data in
terra are SpatVector objects; different from
sf objects.
#load the park boundary
gbr_boundary <- vect("../data/gbr_MPA_boundary.gpkg")
We can plot vector boundaries on top of raster data using
lines():
#plot the SST data and the boundary on top
plot(sst)
lines(gbr_boundary) #this plots the vector lines over the raster data
We now see the outline of the Great Barrier Reef marine park on top of the sea surface temperature data we plotted before.
The data we have at the moment is for a much larger area than just the Great Barrier Reef marine park. To get only that data we need to crop and mask the raster data using the marine park boundary.
Cropping means that we keep only the data inside the extent of the vector we are using. Mask means that all the data outside the vector is set to NA or some other value we specify. Lets have a look how this works.
First lets have a look at the extent of the marine park boundary. We
can get the extent of a raster or vector using ext(). We
need to convert this into a SpatVector object for plotting
using vect(). We only need to do this for plotting; when we
crop, we can just use the marine park boundary as the input.
raster users note: ext() replaces
extent()
gbr_boundary_extent <- ext(gbr_boundary) |>
vect()
plot(sst)
lines(gbr_boundary)
lines(gbr_boundary_extent, col = "blue")
Cropping means we remove everything outside the extent (blue box) of our polygon. Masking sets all values outside our polygon to NA.
So when we crop, we get only the area within the blue box.
We crop using the crop() function, using the raster we
want to crop as the first argument and the vector we are cropping with
second.
#crop
sst_cropped <- crop(sst, gbr_boundary)
#plot
plot(sst_cropped)
lines(gbr_boundary)
Now we have cropped our raster, we can mask it so that we only have
values for the area within the marine park boundary. We do this using
mask:
#mask
sst_cropped_masked <- mask(sst_cropped, gbr_boundary)
#plot
plot(sst_cropped_masked)
lines(gbr_boundary)
Now we only see raster values for cells that are within the marine
park boundary. But remember that the areas that are white, still have
values, they are just NA values.
Often we want to crop and mask one after
the other, and you can do this in one command using
crop(sst, gbr_boundary, mask = TRUE).
For reference, here is a figure comparing what crop,
mask and crop(mask = TRUE) do:
par(mfrow = c(2,2))
plot(sst, main = "Original raster")
lines(gbr_boundary)
plot(sst_cropped, main = "Cropped")
lines(gbr_boundary)
sst |>
mask(gbr_boundary) |>
plot(main = "Masked")
lines(gbr_boundary)
plot(crop(sst, gbr_boundary, mask = TRUE), main = "Cropped and masked")
lines(gbr_boundary)
par(mfrow = c(1,1))
Why not just mask rather than crop and mask? As we see in the figure
above, this would mean we have a lot of area we are not interested in
and even though most of those cells would be NA they take
up space in our raster, so it is not efficient.
Now we have cropped and masked our original raster to get only data
within the area we are interested in, we can start exploring the values.
terra has several functions that can help us do this
easily.
We can get a histogram of all the values in our raster using the
hist() function, which is equivalent to the base R
function.
hist(sst_cropped_masked)
The x-axis is in the units of our raster values; temperature in degrees Celsius in our case. The y-axis is frequency; how many cells in our raster have those values.
You can change the histogram using the same arguments you use with
hist() in base R. For example, lets increase the number of
bars we have:
hist(sst_cropped_masked, breaks = 100)
We can get a frequency table of values in our raster using
freq(). This is essentially the same information that is
shown in the histogram in graphical format.
freq(sst_cropped_masked)
#> layer value count
#> 1 1 20 31
#> 2 1 21 223
#> 3 1 22 396
#> 4 1 23 740
#> 5 1 24 2861
#> 6 1 25 2827
#> 7 1 26 2194
#> 8 1 27 2638
#> 9 1 28 444
The default is to round values to the nearest integer. To get more integers we can do:
freq(sst_cropped_masked, digits = 1) |>
head() #this is a long table: just show the first few values
#> layer value count
#> 1 1 19.8 1
#> 2 1 19.9 3
#> 3 1 20.0 6
#> 4 1 20.1 3
#> 5 1 20.2 3
#> 6 1 20.3 7
We can also check how many NA values are in our raster:
freq(sst_cropped_masked, value = NA)
#> layer value count
#> 1 1 NA 50850
We can use the same summary() command that is used in
base R to get a summary of the statistical information for an entire
raster.
summary(sst_cropped_masked)
#> CRW_SST_2
#> Min. :19.82
#> 1st Qu.:24.25
#> Median :25.05
#> Mean :25.19
#> 3rd Qu.:26.49
#> Max. :27.74
#> NA's :50850
These values are, for example, the mean value of all raster values
excluding NAs.
To get individual statistical values, we need to use
global(). For example to get the mean value of a
raster:
global(sst_cropped_masked, "mean")
#> mean
#> CRW_SST_2 NaN
Hmmm, what went wrong? The default is for global() to
include all values in the raster, and since we have lots of
NAs, the result is NaN. We need to set
na.rm = TRUE to exclude the NA values.
global(sst_cropped_masked, "mean", na.rm = TRUE)
#> mean
#> CRW_SST_2 25.18745
The object returned by global() is a data frame, so if
you want just the value, you need to do:
sst_mean <- as.numeric(global(sst_cropped_masked, "mean", na.rm = TRUE))
sst_mean
#> [1] 25.18745
raster users note: global() replaces
cellStats(), and the default is na.rm = FALSE,
whereas the default for cellStats() was
na.rm = TRUE.
We might want to break our raster values into groups. For example, we
could say that all temperatures that are below the mean temperature are
classified as “cooler” and all temperatures above the mean are “warmer”.
We can do this using the classify function.
raster users note: classify() replaces
reclassify()
#first we create a matrix that will be used for the classification
# all values >= 0 and less than the mean become 1
# all values greater than the mean become 2
reclass_matrix <- c(0, sst_mean, 1,
sst_mean, Inf, 2) |>
matrix(ncol = 3, byrow = TRUE)
#now we classify our raster using this matrix
sst_reclassed <- classify(sst_cropped_masked, reclass_matrix)
#plot the result
plot(sst_reclassed)
A gut check tells us this looks right; the warmer areas are in the north, nearer to the equator.
This plot is ok, but it would be better if the colours were more appropriate and the legend gave some useful information.
plot(sst_reclassed, col = c("blue", "red"), plg = list(legend = c("Cooler", "Warmer")), las = 1) #the las = 1 argument just rotates the y-axis labels so that they are horizontal
The great thing about rasters are you can do maths with them! For
example, doing sst_reclassed + 1 just adds one to each
raster value, and doing sst_reclassed*2 multiplies each
raster value by two.
As an example, lets convert our temperature data into Fahrenheit for our confused colleagues in the U.S. The conversion from Celsius to Fahrenheit is: Fahrenheit = (Celsius * 1.8) + 32.
#do the conversion
sst_fahrenheit <- (sst_cropped_masked*1.8) + 32
#plot our new raster
plot(sst_fahrenheit)
Lets also make the map a bit prettier, choosing more appropriate
colours and adding a scale bar. There are huge number of options for
plotting, see the help file ?help for details. There are
many great plotting packages such as tmap and
ggplot that can be used to make maps, but I will not cover
those here. The Geocomputation
with R website is an excellent resource on map making and geospatial
data in R in general.
plot(sst_fahrenheit, col = hcl.colors(50, palette = "RdYlBu", rev = TRUE), las = 1) #see ?hcl.colors for more info on colour palettes. rev= TRUE because we want to reverse the palette: red colours are the highest values and blue the lowest
lines(gbr_boundary)
sbar(d = 400, type = "bar", divs = 2, below = "km") #400km scale bar with 2 divisions and "km" written below
The swatches below show all the colour scales available via
hcl.colors(). The code is directly from
?hcl.colors().
## color swatches for HCL palettes
hcl.swatch <- function(type = NULL, n = 5, nrow = 11,
border = if (n < 15) "black" else NA) {
palette <- hcl.pals(type)
cols <- sapply(palette, hcl.colors, n = n)
ncol <- ncol(cols)
nswatch <- min(ncol, nrow)
par(mar = rep(0.1, 4),
mfrow = c(1, min(5, ceiling(ncol/nrow))),
pin = c(1, 0.5 * nswatch),
cex = 0.7)
while (length(palette)) {
subset <- 1:min(nrow, ncol(cols))
plot.new()
plot.window(c(0, n), c(0, nrow + 1))
text(0, rev(subset) + 0.1, palette[subset], adj = c(0, 0))
y <- rep(subset, each = n)
rect(rep(0:(n-1), n), rev(y), rep(1:n, n), rev(y) - 0.5,
col = cols[, subset], border = border)
palette <- palette[-subset]
cols <- cols[, -subset, drop = FALSE]
}
par(mfrow = c(1, 1), mar = c(5.1, 4.1, 4.1, 2.1), cex = 1)
}
hcl.swatch("qualitative")
hcl.swatch("sequential")
hcl.swatch("diverging")
hcl.swatch("divergingx")
We have been looking at data within the Great Barrier Reef Marine Park boundary. What if we want to extract data for some areas within that area? For example we might want to know the temperature within some of the zones within the marine park. We could crop and mask our data for each of our zones, but this could get messy if we have a lot of zones. Happily there is an easier way.
First lets load the geospatial data for two habitat protection zones in the Great Barrier Reef. This is just a small part of the Great Barrier Reef Marine Park zoning map. Surprisingly, habitat protection zones allow all types of fishing except trawling. You can view more details about the zones and the activities that are allowed here.
#load the zones data
zones <- vect("../data/gbr_habitat_protection_zones.gpkg")
#take a look
zones
#> class : SpatVector
#> geometry : polygons
#> dimensions : 2, 19 (geometries, attributes)
#> extent : 145.9373, 152.116, -22.41511, -16.38992 (xmin, xmax, ymin, ymax)
#> source : gbr_habitat_protection_zones.gpkg
#> coord. ref. : lon/lat GDA2020 (EPSG:7844)
#> names : LEG_NAME DATE SCHED_NO SCHED_DESC PART_NO
#> type : <chr> <chr> <int> <chr> <int>
#> values : Great Barrier ~ Gazetted 2004 1 Amalgamated Gr~ 2
#> Great Barrier ~ Gazetted 2004 1 Amalgamated Gr~ 2
#> TYPE NAME GROUP_ID CODE PERMIT_DESC (and 9 more)
#> <chr> <chr> <chr> <int> <chr>
#> Habitat Protec~ HP-16-5126 HPZ 10 Arlington Reef~
#> Habitat Protec~ HP-21-5296 HPZ 10 Capricorn Channel
There are lots of columns, but lets use the “PERMIT_DESC” column, which is the zone location, to map them since it has distinct and interesting descriptions for them.
#plot the zones with our temperature data
plot(zones, "PERMIT_DESC")
Lets plot these zones over our temperature data for context.
plot(sst_cropped_masked)
lines(zones)
Geographic coordinate systems: uses a 3-D surface to define locations on the Earth using longitude and latitude
Projected coordinate reference system: translates a GCS into 2-D so we can measure distances, areas, etc. in more useful units such as meters or, dare I say, feet.
https://gis.stackexchange.com/questions/149749/is-wgs84-a-coordinate-system-or-projection-system
https://r.geocompx.org/spatial-class#crs-intro
https://docs.qgis.org/3.34/en/docs/gentle_gis_introduction/coordinate_reference_systems.html
Although we have been able to plot the zones and the temperature
raster together, they are actually in different coordinate reference
systems (crs’s). To see what the crs of a raster or vector is we use
crs().
crs(sst_cropped_masked)
#> [1] "GEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]]"
Well, that’s a lot of information. This is the “well-known text” version of the coordinate reference system. To get just some basic information that most of the time is all we need, we can do:
crs(sst_cropped_masked, describe = TRUE)
#> name authority code area extent
#> 1 WGS 84 EPSG 4326 <NA> NA, NA, NA, NA
crs(zones, describe = TRUE)
#> name authority code
#> 1 GDA2020 EPSG 7844
#> area
#> 1 Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore
#> extent
#> 1 93.41, 173.34, -8.47, -60.55
zonal(sst_cropped_masked, zones, fun = "mean")
#> CRW_SST_2
#> 1 26.11176
#> 2 24.26957
zonal(sst_cropped_masked, zones, fun = "max")
#> CRW_SST_2
#> 1 26.71
#> 2 24.68
zones_data_df <- extract(sst_cropped_masked, zones)
#> Warning: [extract] transforming vector data to the CRS of the raster
zones_data_df
#> ID CRW_SST_2
#> 1 1 26.71
#> 2 1 26.61
#> 3 1 26.63
#> 4 1 26.51
#> 5 1 25.95
#> 6 1 26.16
#> 7 1 26.39
#> 8 1 26.60
#> 9 1 25.59
#> 10 1 25.79
#> 11 1 26.00
#> 12 1 26.25
#> 13 1 25.42
#> 14 1 25.63
#> 15 1 25.85
#> 16 1 26.12
#> 17 1 25.69
#> 18 2 23.97
#> 19 2 24.04
#> 20 2 24.11
#> 21 2 24.17
#> 22 2 24.21
#> 23 2 24.24
#> 24 2 24.29
#> 25 2 24.34
#> 26 2 24.39
#> 27 2 24.44
#> 28 2 23.92
#> 29 2 24.00
#> 30 2 24.06
#> 31 2 24.10
#> 32 2 24.14
#> 33 2 24.16
#> 34 2 24.21
#> 35 2 24.25
#> 36 2 24.30
#> 37 2 24.35
#> 38 2 24.40
#> 39 2 24.45
#> 40 2 23.92
#> 41 2 23.98
#> 42 2 24.02
#> 43 2 24.05
#> 44 2 24.07
#> 45 2 24.08
#> 46 2 24.12
#> 47 2 24.15
#> 48 2 24.20
#> 49 2 24.26
#> 50 2 24.32
#> 51 2 24.38
#> 52 2 24.44
#> 53 2 23.94
#> 54 2 23.97
#> 55 2 24.00
#> 56 2 24.01
#> 57 2 24.01
#> 58 2 24.01
#> 59 2 24.04
#> 60 2 24.07
#> 61 2 24.12
#> 62 2 24.18
#> 63 2 24.25
#> 64 2 24.32
#> 65 2 24.39
#> 66 2 24.45
#> 67 2 24.50
#> 68 2 23.94
#> 69 2 23.97
#> 70 2 23.98
#> 71 2 23.98
#> 72 2 23.97
#> 73 2 23.97
#> 74 2 23.99
#> 75 2 24.03
#> 76 2 24.08
#> 77 2 24.15
#> 78 2 24.22
#> 79 2 24.30
#> 80 2 24.38
#> 81 2 24.44
#> 82 2 24.48
#> 83 2 24.52
#> 84 2 24.56
#> 85 2 23.93
#> 86 2 23.97
#> 87 2 23.99
#> 88 2 23.99
#> 89 2 23.98
#> 90 2 23.99
#> 91 2 24.01
#> 92 2 24.05
#> 93 2 24.10
#> 94 2 24.18
#> 95 2 24.26
#> 96 2 24.34
#> 97 2 24.42
#> 98 2 24.47
#> 99 2 24.50
#> 100 2 24.53
#> 101 2 24.56
#> 102 2 24.58
#> 103 2 24.60
#> 104 2 23.88
#> 105 2 23.95
#> 106 2 23.99
#> 107 2 24.02
#> 108 2 24.02
#> 109 2 24.05
#> 110 2 24.08
#> 111 2 24.12
#> 112 2 24.17
#> 113 2 24.25
#> 114 2 24.34
#> 115 2 24.42
#> 116 2 24.49
#> 117 2 24.53
#> 118 2 24.55
#> 119 2 24.57
#> 120 2 24.59
#> 121 2 24.61
#> 122 2 24.62
#> 123 2 24.64
#> 124 2 23.82
#> 125 2 23.92
#> 126 2 24.00
#> 127 2 24.05
#> 128 2 24.08
#> 129 2 24.12
#> 130 2 24.16
#> 131 2 24.20
#> 132 2 24.26
#> 133 2 24.33
#> 134 2 24.42
#> 135 2 24.50
#> 136 2 24.57
#> 137 2 24.60
#> 138 2 24.61
#> 139 2 24.62
#> 140 2 24.63
#> 141 2 24.64
#> 142 2 24.65
#> 143 2 24.66
#> 144 2 24.66
#> 145 2 24.64
#> 146 2 23.77
#> 147 2 23.89
#> 148 2 24.00
#> 149 2 24.08
#> 150 2 24.13
#> 151 2 24.18
#> 152 2 24.22
#> 153 2 24.25
#> 154 2 24.31
#> 155 2 24.36
#> 156 2 24.44
#> 157 2 24.52
#> 158 2 24.58
#> 159 2 24.61
#> 160 2 24.63
#> 161 2 24.64
#> 162 2 24.65
#> 163 2 24.66
#> 164 2 24.67
#> 165 2 24.67
#> 166 2 23.72
#> 167 2 23.85
#> 168 2 23.98
#> 169 2 24.07
#> 170 2 24.13
#> 171 2 24.19
#> 172 2 24.22
#> 173 2 24.26
#> 174 2 24.30
#> 175 2 24.35
#> 176 2 24.42
#> 177 2 24.49
#> 178 2 24.54
#> 179 2 24.58
#> 180 2 24.60
#> 181 2 24.63
#> 182 2 24.65
#> 183 2 24.67
#> 184 2 24.68
#> 185 2 24.68
#> 186 2 23.67
#> 187 2 23.81
#> 188 2 23.94
#> 189 2 24.04
#> 190 2 24.10
#> 191 2 24.14
#> 192 2 24.17
#> 193 2 24.21
#> 194 2 24.25
#> 195 2 24.30
#> 196 2 24.36
#> 197 2 24.42
#> 198 2 24.47
#> 199 2 24.51
#> 200 2 24.55
#> 201 2 24.59
#> 202 2 24.63
#> 203 2 24.66
#> 204 2 24.68
#> 205 2 24.68
boxplot(CRW_SST_2 ~ ID, zones_data_df)
A very useful feature of rasters is that they can have many layers.
These layers often represent different time periods, such as days or
months. Lets look at a multi-layer raster which has the same temperature
data that we have been looking at, but each layer represents the mean
temperature in a month. We load the multi-layer raster in the same way
as any other raster, using rast():
sst_monthly <- rast("../data/gbr_monthly_temp.tif")
sst_monthly
#> class : SpatRaster
#> dimensions : 301, 261, 36 (nrow, ncol, nlyr)
#> resolution : 0.05, 0.05 (x, y)
#> extent : 142, 155.05, -25, -9.95 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)
#> source : gbr_monthly_temp.tif
#> names : ym_202006, ym_202007, ym_202008, ym_202009, ym_202010, ym_202011, ...
#> min values : 19.80000, 18.51452, 19.00613, 21.096, 23.04516, 24.715, ...
#> max values : 29.47767, 28.65097, 28.51323, 28.864, 29.41645, 29.452, ...
We can that the raster has 36 layers. Lets plot it to see what we get.
plot(sst_monthly)
We get one map for each raster layer, but not all the data because
you wouldn’t be able to see the maps if they were much smaller. We can
plot a specific subset of the data using double square brackets
[[]] to select a set of layers. For example, if we wanted
the last 4 layers in the raster, which are layers 33 to 36:
plot(sst_monthly[[33:36]])
Happily, we can use the same functions we have learnt on this multi-layer raster. For example, we can convert all the rasters to Fahrenheit:
sst_monthly_fahrenheit <- (sst_monthly * 1.8) + 32
plot(sst_monthly_fahrenheit)
We can also crop and mask our data:
sst_monthly_cropped_masked <- crop(sst_monthly, gbr_boundary, mask = TRUE)
plot(sst_monthly_cropped_masked)